The CSEM technique is an important geophysical tool for hydrocarbon prospecting in the earth's subsurface. In a CSEM survey, an electromagnetic-wave source (transmitter) generates an electromagnetic wave. The electromagnetic signal induced in the earth by the transmitter is recorded constantly in time by one or more receivers. The electromagnetic signal at a receiver location depends on physical properties, especially the electrical properties, of the medium in which the electromagnetic wave has passed through from the source to the receiver. The behavior of this signal as a function of frequency and transmitter location or separation (offset) between transmitter and receiver can be used to estimate the spatially varying resistivity model of the subsurface within a certain depth range. This estimated subsurface resistivity model is used for identifying resistivity anomalies indicating the presence of hydrocarbons (oil or gas) in the earth's subsurface.
FIG. 1 illustrates a typical marine CSEM survey in which a constantly active electromagnetic-wave transmitter 11 is towed below the water surface 15 along a line 12 above electromagnetic receivers 13 (two neighboring receivers are shown) deployed on the seafloor 14. Reference number 16 indicates the offset between the right-most receiver and the source when the source is at location 11A. For more details see Chapter 12, page 931 in Investigations In Geophysics No. 3, Electromagnetic Methods In Applied Geophysics, volume 2, edited by Misac N. Nabighian, Society of Exploration Geophysicists, 1991). Alternative configurations include stationary transmitters on the seafloor or in the water column as well as magnetic transmitter antennae and connecting several receivers in a towed array (see, for example, U.S. Pat. No. 4,617,518 to Srnka). The receivers typically have multiple sensors designed to record different vector components of the electric and/or magnetic fields. A sensor is also called a channel. The data recorded in one channel correspond to one vector component of the electromagnetic field. Every receiver records the electromagnetic signal constantly in time during a survey. The data recorded by one sensor at a receiver location are normally called a common-receiver gather, or simply called a receiver gather. Under the stationary-receiver configuration, a common-receiver gather represents the electromagnetic signal at the fixed receiver location induced by the source at all different source locations, or at different times during the survey. Similarly, data can also be sorted in common-source gathers to represent the electromagnetic field at those receiver locations from a source at a fixed source location.
Marine CSEM data are typically interpreted in the temporal frequency domain. After taking out the frequency-dependent effects of the source and the receiver themselves, the signal at a frequency represents the response of the earth to electromagnetic signal at that temporal frequency. It is this response that provides us information about the subsurface electrical properties. Like any other type of wave, the electromagnetic signal in a CSEM survey has two attributes, amplitude and phase. The signals are therefore conveniently represented as complex numbers in either rectangular (real-imaginary) or polar (amplitude-phase) form.
In practice, the receiver data are usually converted to temporal frequency by dividing (or “binning”) the recorded time-domain data into time intervals (i.e. bins: x1, x2, and x3 as shown in FIG. 2A) and determining the spectrum within each bin by standard methods based on the Fourier Transform. The signal recorded by a receiver is 21 and, for reference, the transmitted periodic waveform 22 is also shown. FIG. 2B shows the amplitudes of the spectral components from the bin x3. Unlike the example shown in FIG. 2A, a typical bin length is several periods of the transmitter waveform. Each bin might correspond to a different position of the source arrow 11 in FIG. 1. Some methods of transforming data to the time-frequency domain include the Short-Time Fourier Transform (J. Allen, L. Rabiner, “A Unified Approach to Short-Time Fourier Analysis and Synthesis,” Proc. of the IEEE 65, 1558-64, (1977)); and the Choi-Williams transform (H. Choi and W. Williams, “Improved time-frequency representation of multicomponent signals using exponential kernels,” IEEE Trans. on Acoust., Speech, and Signal Processing, 37, 862-871, (1989)). In the temporal-frequency domain, signals recorded by a receiver, including both amplitude and phase, of each of the temporal-frequency components are functions of bin, or the transmitter location, or the signed offset distance between source and receiver. FIGS. 3A-B show an example of amplitude (3A) and phase (3B) variation versus transmitter-receiver offset at frequency ⅜ Hz. The drawings represent model calculations with the solid-line curves representing a resistivity model containing a hydrocarbon (high resistivity) layer, whereas the dashed line curves were generated using a resistivity model without a reservoir. As shown in the drawings, both the phase and amplitude of CSEM data can be indicative of resistive (and potentially hydrocarbon-bearing) strata in the subsurface, and can thus be used to estimate the subsurface rock electrical conductivity or resistivity. Hydrocarbon bearing rocks usually show higher resistivity than the surrounding sediments. The differences between the solid-line curves and the dashed-line curves show how CSEM data may be used to detect the presence of hydrocarbons. Thus, the subsurface rock resistivity information derived from the CSEM data is valuable for hydrocarbon exploration risk reduction.
The estimation of the subsurface resistivity (or conductivity) model in three-dimensional (3-D) space from measured CSEM data is an inverse problem. Solving an inverse problem is a trial-and-error iterative process. The final estimated model should be able to predict data that match the measured data and satisfy any constraints that may be applicable to the model.
This process (i.e. updating the resistivity model for the next iteration) can be either human-guided manual adjustment of the subsurface resistivity model or an automatic model update predicted from some appropriate mathematical measures of the misfit between measured and the predicted data. See for example, G. A. Newman and D. L. Alumbaugh, “Three-dimensional massively parallel electromagnetic inversion—I. Theory,” Geophys. J. Int., 128, 345-354 (1997) and Y. Sasaki, “Full 3-D inversion of electromagnetic data on PC,” J. of Applied Geophys., 46, 45-54, (2001), or a combination of the two. The prediction of electromagnetic data from a resistivity model of the subsurface is achieved by numerically solving Maxwell's electromagnetic field equations, a process called forward modeling.
In many examples of CSEM hardware, data cannot be effectively recorded at the nearest offsets because the dynamic range of the receiver's digitizers is too small to accommodate the large dynamic range of the data. This region is sometimes known as the “saturation zone” and typically encompasses source-receiver offsets of less than 500 meters depending on amplifier property of the receiver. An example is shown in FIG. 3A, in which constant amplitude is observed for offset roughly within 500 m.
The inversion of CSEM data for the subsurface conductivity is a computationally intensive process, since it involves many forward simulations of the electromagnetic field in multi-dimensional space. To speed up the inversion process in multi-dimensional space, such as 2-D or 3D space, the model-update prediction is derived from the forward modeling and the transmitter-receiver reciprocity property can be used to reduce the number of forward modeling operations; see the previously cited Newman and Alumbaugh reference, and also U.S. Provisional Patent Application No. 60/780,232. By using the reciprocity principle (switching the role of a transmitter and a receiver), the electromagnetic fields in one entire receiver gather (as shown in FIGS. 3A-B) can be obtained in one forward modeling operation by calculating the electromagnetic fields at those original transmitter locations from a transmitter located at the original receiver location. In the traditional frequency-domain inversion process, every receiver gather needs to be forward simulated separately and compared to the measured data. The number of receiver gathers to be simulated in a survey is the product of the number of receivers deployed, frequencies to be used, and the number of components of each receiver. If the inversion is performed in time domain, the number of forward simulations is proportional to the number of receiver or source gathers in a survey depending on whether the transmitter-receiver reciprocity is used. The large number of independent forward simulations required in the current inversion techniques discourages its application to large 3-D surveys, such as reconnaissance surveys covering large area with regularly spaced receivers as illustrated by an example in FIG. 4. Techniques leading to a substantial speedup of CSEM data inversion in multi-dimensional space are crucial for its application in 3-D surveys. The present invention fulfills this need.
FIG. 4 shows a surface map view of a typical data acquisition pattern for a CSEM reconnaissance survey. The distribution of receivers (indicated by both black and white arrows) on a regular 2-D grid is shown as well as a set of parallel transmitter towlines (dashed lines); the receiver interval is usually several kilometers. Acquiring data along a regular grid is natural for CSEM reconnaissance, where a priori information about the subsurface is limited.
Current methods of CSEM data inversion will next be examined in somewhat more detail. As previously stated, CSEM data inversion is an iterative method for determining the resistivity of the subsurface from CSEM data measured at the earth's surface or seafloor. The result of inversion is a geo-electric model of the subsurface obtained by updating a starting model of the earth resistivity to minimize the mismatch between measured and simulated data. The model update from iteration to iteration can be achieved by either human-guided manual adjustment of the resistivity model or an automatic model update predicted from some appropriate mathematical measures of the misfit between measured and the predicted data (see for example, G. A. Newman and D. L. Alumbaugh, op. cit.) or a combination of the two.
Most of the geological and electrical information that may be available about the subsurface, such as structural and rock physical property information from seismic data and electrical property information from available well measurements, can be taken into account by human-guided manual model updates more easily than by automatic model updates predicted from some appropriate mathematical measures of the mismatch between measured and the predicted data. However, the human-guided manual model update becomes awkward as the survey size and/or the subsurface geology complexity increases. This is especially true for 3D inversion due to the great flexibility in updating the model in 3-D space.
Most of the current inversion procedures adopt some automatic model-update schemes based on numerical optimization procedures which adjust subsurface resistivities and possibly other parameters until the defined objective function is reduced to a sufficiently small value. The objective function usually includes term(s) describing the data mismatch between the forward simulated data and the measured data and other term(s) describing some geological information inputs and constraints. Some model constraints may also be enforced directly in the model-update process. The inversion process drives the model along the direction of reducing the data mismatch and satisfying any geological constraints included either in the objective function or enforced in the model-update process.
These inverted models from either manual or automatic model-update processes should be able to produce synthetic CSEM data that accurately match the measured data. Inversions using either of the model-update procedures outlined above require repeated solution of Maxwell's equations (or forward modeling) for a large number of models and transmitter-receiver configurations. The forward modeling of CSEM data in 3D space is computationally intensive and it dominates the computational time and costs in the CSEM data inversion (see, for example, D. L. Alumbaugh and G. A. Newman, “3-D massively parallel electromagnetic inversion—Part II, Analysis of a cross well experiment,” Geophys. J. Int. 128, 355-363 (1997)). Under some simple situations, the subsurface resistivities might be approximated by a 1-D layered model which limits any variation of resistivity along the horizontal direction for more efficient forward modeling and inversion (see, for example, S. Constable and C. J. Weiss, “Mapping thin resistors and hydrocarbons with marine EM methods: Insights from ID modeling”, Geophysics 71, G43-G51, (2006)). In general, such simplification is not accurate enough for application in hydrocarbon exploration.
There exist several forward-modeling schemes for the simulation of electromagnetic wave propagation. The commonly-used ones for general 3-D models are the finite-difference method, hereinafter “FDM”, the finite element method, hereinafter “FEM”, and the integral equation method, hereinafter “IEM.” These are the standard approaches for numerically solving any partial differential equation (s) that cannot be solved analytically. In practical applications of these methods, the physical properties, such as the resistivity and dielectric permittivity, are represented by discrete cells in the whole space of interest, or in a localized zone for some special applications of the IEM. The governing equations, Maxwell's equations for CSEM applications, are represented in discrete forms on the cell grids for both FDM and FEM and are used to solve the electromagnetic field numerically over the cell grids. The FDM normally uses rectangular cells without assuming any particular geometric structures of the physical property in space (G. A. Newman and D. L. Alumbaugh, op. cit.). The FEM normally uses more general geometric shapes than rectangles (J. H. Coggon, “Electromagnetic and Electrical Modeling by Finite Element Method”, Geophysics 36, 132-155 (1971)) that are able to represent the model in more detail than the FDM at the expense of more complex model representation and governing equations over the cell grids.
The IEM recasts the system of differential equations implied by Maxwell's equations into an associated integral equation by making use of the properties of the Green's function for the electric and/or magnetic field in a uniform or layered model. A uniform or layered material is typically used for the reference Green's function because highly accurate and rapidly computed solutions are available for these models. The resulting integral equations naturally give rise to computational schemes that work very well for compact objects imbedded in a uniform or layered background (such as a ship in the deep ocean or an aircraft high in the atmosphere).
The forward modeling methods described in the preceding paragraphs can be applied in both time and frequency domain (for time domain example, see for example, M. Commer and G. A. Newman, “A parallel finite-difference approach for 3D transient electromagnetic modeling galvanic sources”, Geophysics. 69, 1192-1202 (2004)). Forward modeling CSEM data in the time domain offers advantages in handling the so-called air-wave effect in land or shallow water surveys (the air wave is the direct transmission from the broadcasting antenna to the detecting antenna through the air). However, it is computationally more costly than in frequency domain due to the large number of time steps needed to simulate the propagation of electromagnetic waves in the model.
All of the preceding forward-modeling approaches in frequency domain result in a very large linear system to solve. The large size of the linear system combined with the large number of forward modelings needed for a survey makes the forward modeling time consuming. A powerful computer is often needed in order to obtain results in a reasonable time. A number of techniques have been developed to speed up the computation at different stages of the inverse process. For example, more efficient optimization techniques such as the non-linear conjugate gradient (NLCG) solver, multi-grid solvers, approximate computation of the sensitivity matrix, source-receiver reciprocity, etc. All those techniques are helpful, but more improvements are needed to make electromagnetic inversion in 3D space a routine practice with reasonable computer resources.